\documentclass[a4paper]{article} 

\usepackage{listings}
\lstset{tabsize=4, %
	frame=shadowbox, %把代码用带有阴影的框圈起来
	commentstyle=\color{red!50!green!50!blue!50},%浅灰色的注释
	rulesepcolor=\color{red!20!green!20!blue!20},%代码块边框为淡青色
	keywordstyle=\color{blue!90}\bfseries, %代码关键字的颜色为蓝色，粗体
	showstringspaces=false,%不显示代码字符串中间的空格标记
	stringstyle=\ttfamily, % 代码字符串的特殊格式
	keepspaces=true, %
	breakindent=22pt, %
	numbers=left,%左侧显示行号
	stepnumber=1,%
	numberstyle=\tiny, %行号字体用小号
	basicstyle=\footnotesize, %
	showspaces=false, %
	flexiblecolumns=true, %
	breaklines=true, %对过长的代码自动换行
	breakautoindent=true,%
	breakindent=4em, %
	%escapebegin=\begin{CJK*}{GBK}{hei},escapeend=\end{CJK*},
	aboveskip=1em, %代码块边框
	fontadjust,
	captionpos=t,
	framextopmargin=2pt,framexbottommargin=2pt,abovecaptionskip=-3pt,belowcaptionskip=3pt,
	xleftmargin=4em,xrightmargin=4em, % 设定listing左右的空白
	texcl=true,
	% 设定中文冲突，断行，列模式，数学环境输入，listing数字的样式
	extendedchars=false,columns=flexible,mathescape=true
	% numbersep=-1em
}
\usepackage{cite}
\usepackage{fancyhdr}
\usepackage{lipsum}
\usepackage{extarrows}
\usepackage{ctex}
\usepackage{amsfonts,amssymb} 
\usepackage{amsmath} 
\usepackage{cases}
\usepackage{geometry}
\usepackage{graphicx}
\usepackage{listings}
\usepackage{xcolor}
\usepackage{framed}
\usepackage{float}
\usepackage{ulem}
\usepackage{booktabs}
\usepackage{amsthm}
\usepackage{tikz}
\usepackage{pgfplots}
\pgfplotsset{width=10cm,compat=1.9}
\geometry{a4paper,left=3cm,right=3cm,top=3cm,bottom=3cm}

\begin{document}
	\begin{center}
		{\heiti {\huge 第二次编程作业报告}}
	\end{center}
	
	\begin{center}
		{\large 求数2101\ 叶陈昊\ 3210106359}
	\end{center}

	\section{ProblemA}
	本题要求定义一个函数，将插值点组$x_0, x_1, \cdots, x_n$、函数$f$以及实数$x$作为参数传入，将函数$f$的插值多项式$p$在$x$处的值返回.
	
	可以发现，这个函数仅在题B、C中有用处，而在题D、E中是没有用处的（因为题D、E中并没有给出明确的函数$f$，而是给出$f(x_0), f(x_1)\cdots, f(x_n)$这样的离散值）.因此，为了使这个牛顿插值多项式更具有适用性，我们另写一个头文件\texttt{NewtonPoly.h}，来定义一些比较具有适用性的函数.
	
	对于问题A，我们在\texttt{NewtonPoly.h}中声明如下函数：
	\begin{lstlisting}[language=c++]
		double Newton_poly_at_x(const vector<double> &inter_pts, const vector<double> &vals_of_pts, double x);
	\end{lstlisting}
	
	这个函数传入了插值点组，对应的函数值组以及一个实数，返回的是由这些点构成的插值多项式在这个实数处的值.
	
	根据这个函数，我们很容易想到题A的处理办法，即将$f(x_0), f(x_1), \cdots, f(x_n)$存到一个数组里，作为传入\texttt{Newton\_poly\_at\_x}的函数值组参数.具体见\texttt{ProblemA.h}、\texttt{ProblemA.cpp}.
	
	另外，\texttt{NewtonPoly.h}中还声明了一些函数，在做后面题的时候会有所帮助：
	
	\begin{itemize}
		\item \texttt{Newton\_poly} ：得到牛顿插值多项式的表达式，可以用于在软件上输入多项式直接作图（但多项式次数太高亲测作图会产生一些精度问题，最后还是选择了描点连线法.不过写都写了，就放在这里叭）.
		\item \texttt{Newton\_poly\_at\_x}：得到牛顿插值多项式在$x$处的值，具体上面已讲.
		\item \texttt{Newton\_dif\_at\_x}：得到牛顿插值多项式在$x$处的一次导数的值，具体细节我们在题D讨论.
		\item \texttt{get\_diff\_table}：得到差商表（含广义差商表，详细见题D）.
		\item \texttt{print\_diff\_table}：打印差商表.
	\end{itemize}
	
	\section{ProblemB}
	
	问题B给定函数$f(x)=\frac{1}{1+x^2}$，要求我们在$[-5,5]$上取$n+1$个等距节点并依据这些点绘制$n$次插值多项式的图像.我们使用\texttt{ProblemA.h}中声明的函数\texttt{poly\_at\_x}，进行“描点连线”.等距取201个采样点，对每个采样点$x$，将函数$f$、插值点$x_i,i=0,1,\cdots,n$，和$x$传入\texttt{poly\_at\_x}得到插值多项式在$x$处的值$y$，之后描点$(x,y)$即可.另外，将这201个点的横纵坐标输出到\texttt{txt}文件中，方便后面使用\LaTeX 作图.
	
	对于$n=2,4,6,8$，如此操作之后，我们在区间$[-5,5]$上作出这些插值多项式的图像，以及函数$f$的图像，如下图所示.可以发现，当$n$越来越大时，离原点较远处，多项式的振荡越来越明显，因此出现Runge现象.
	
	%第一次尝试直接传多项式，在这里多项式次数不高，作图还是较精确的
	\begin{center}
		\begin{tikzpicture}
			\begin{axis}[
				title={题B：等距节点产生Runge现象},
				legend pos=outer north east,
				xlabel=$x$,
				ylabel={$y$},
				]
				
				%作图f(x)-绿色
				\addplot [
				domain=-5:5, 
				samples=100, 
				color=green,
				]
				{1/(1+x^2)};
				\addlegendentry{$f(x)=\frac{1}{1+x^2}$}
				
				%2次插值-红色
%				\addplot [
%				domain=-5:5, 
%				samples=100, 
%				color=red,
%				]
%				{0.0384615+0.192308*(x+5)-0.0384615*(x+5)*x};
				\addplot[red, mark=none] table {ProblemB_data/file1.txt};
				\addlegendentry{$n=2$}
				
				%4次插值-蓝色			
%				\addplot [
%				domain=-5:5, 
%				samples=100, 
%				color=blue,
%				]
%				{0.0384615+0.0397878*(x+5)+0.061008*(x+5)*(x+2.5)-0.0265252*(x+5)*(x+2.5)*x+0.00530504*(x+5)*(x+2.5)*x*(x-2.5)};
				\addplot[blue, mark=none] table {ProblemB_data/file2.txt};
				\addlegendentry{$n=4$}
				
				%6次插值-黑色				
%				\addplot [
%				domain=-5:5, 
%				samples=100, 
%				color=black,
%				]
%				{0.0384615+0.0264644*(x+5)+0.0248454*(x+5)*(x+3.33333)+0.0149446*(x+5)*(x+3.33333)*(x+1.66667)-0.0131699*(x+5)*(x+3.33333)*(x+1.66667)*x+0.00420316*(x+5)*(x+3.33333)*(x+1.66667)*x*(x-1.66667)-0.000840633*(x+5)*(x+3.33333)*(x+1.66667)*x*(x-1.66667)*(x-3.33333)};
				\addplot[black, mark=none] table {ProblemB_data/file3.txt};
				\addlegendentry{$n=6$}
				
				%8次插值-橙色				
%				\addplot [
%				domain=-5:5, 
%				samples=100, 
%				color=orange,
%				]
%				{0.0384615+0.0223428*(x+5)+0.013956*(x+5)*(x+3.75)+0.0117043*(x+5)*(x+3.75)*(x+2.5)+0.000674338*(x+5)*(x+3.75)*(x+2.5)*(x+1.25)-0.00489646*(x+5)*(x+3.75)*(x+2.5)*(x+1.25)*x+0.00243964*(x+5)*(x+3.75)*(x+2.5)*(x+1.25)*x*(x-1.25)-0.000687223*(x+5)*(x+3.75)*(x+2.5)*(x+1.25)*x*(x-1.25)*(x-2.5)+0.000137445*(x+5)*(x+3.75)*(x+2.5)*(x+1.25)*x*(x-1.25)*(x-2.5)*(x-3.75)};
				\addplot[orange, mark=none] table {ProblemB_data/file4.txt};
				\addlegendentry{$n=8$}
				
			\end{axis}
		\end{tikzpicture}
	\end{center}
	
	\section{ProblemC}
	
	问题C给定函数$f(x)=\frac{1}{1+25x^2}$，要求我们在$[-1,1]$上取$n$个点（Chebyshev多项式$T_n$
	的$n$个零点）并依据这些点绘制$n$次插值多项式的图像.这跟题B的思路是一致的.我们使用\texttt{ProblemA.h}中声明的函数\texttt{poly\_at\_x}，进行“描点连线”.等距取201个采样点，对每个采样点$x$，将函数$f$、插值点$x_i,i=0,1,\cdots,n$，和$x$传入\texttt{poly\_at\_x}得到插值多项式在$x$处的值$y$，之后描点$(x,y)$即可.另外，将这201个点的横纵坐标输出到\texttt{txt}文件中，方便后面使用\LaTeX 作图.
	
	对于$n=5,10,15,20$，如此操作之后，我们在区间$[-1,1]$上作出这些插值多项式的图像，以及函数$f$的图像，如下图所示.可以发现，当$n$越来越大时，插值多项式的拟合效果还是不错的，没有在部分区域上振荡剧烈，基本上都是均匀的.
	
	%尝试直接传多项式，在n=20时，图像较精确图有所误差，故不用这种方法了
	\begin{center}
		\begin{tikzpicture}
			\begin{axis}[
				title=题C：非等距节点的插值可行性,
				legend pos=outer north east,
				xlabel = $x$,
				ylabel = {$y$},
				]
				
				%函数f(x)-绿色
				\addplot [
				domain=-1:1, 
				samples=100, 
				color=green,
				]
				{1/(1+25*x^2)};
				\addlegendentry{$f(x)=\frac{1}{1+25x^2}$}

				%n=5—红色
				\addplot[red, mark=none] table {ProblemC_data/file1.txt};
				\addlegendentry{$n=5$}
				
				%n=10—蓝色
				\addplot[blue, mark=none] table {ProblemC_data/file2.txt};
				\addlegendentry{$n=10$}
				
				%n=15—黑色
				\addplot[black, mark=none] table {ProblemC_data/file3.txt};
				\addlegendentry{$n=15$}
				
				%n=20—橙色 
				\addplot[orange, mark=none] table {ProblemC_data/file4.txt}; 
				\addlegendentry{$n=20$}
				
			\end{axis}
		\end{tikzpicture}
	\end{center}
	
	\section{ProblemD}
	
	问题D相当于给出了插值点$t_0=0,t_1=3,t_2=5,t_3=8,t_4=13$，位移（相应函数值）$s=s(t),s_0=0,s_1=225,s_2=383,s_3=623,s_4=993$，速度（相应导数值）$v=v(t)=s'(t),v_0=75,v_1=77,v_2=80,v_3=74,x_4=72.$
	并依据这些数据绘制插值多项式的图像.由于这题有导数值的存在，是Hermite插值.不过根据课上所讲内容，其有一种解法和Newton插值的解法是差不多的，即作出差商表，只不过Hermite插值需要广义差商作支撑.
	
	而广义差商表的制作不是很容易.为了简便，我们需要用户将插值点按从小到大的方式输入数组，并将含有导数值的插值点重复输入一次（有几个就输入几个）；而函数及其导数值一并输入另一个数组，要求同节点的靠在一起，并按导数阶数依次排列.这样插值点个数和函数（导数）值的个数相等，我们不仅传入了所有数据，而且还将这种特殊情形转化成了一般情形（即插值点互不相同的情形）.我们可以举一些例子让用户了解如何传入这些数组.
	\begin{itemize}
		\item 理论作业第5题：
		
		传入插值点：$\{0,1,1,1,2,2\}$，
		
		传入函数（导数）值：$\{f(0),f(1),f'(1),f''(1)/2,f(2),f'(2)\}=\{0,1,7,42/2=21,128,448\}$.
		\item 理论作业第6题：
		
		传入插值点：$\{0,1,1,3,3\}$，
		
		传入函数（导数）值：$\{f(0),f(1),f'(1),f(3),f'(3)\}=\{1,2,-1,0,0\}$.
	\end{itemize}
	
	我们转化完思路后，接下来就是改进差商表的过程，即改进在\texttt{NewtonPoly.h}中声明的函数\texttt{get\_diff\_table}.这部分见代码即可.
	
	还有一点是，题目要求输出插值多项式在$t=10$时的函数值（位移$s$）和导数值（速度$v$），函数值仍然可以通过（\texttt{cpp}中的）函数\texttt{Newton\_poly\_at\_x}实现，但导数值就没有现有的（\texttt{cpp}中的）函数可用了.我们得在\texttt{NewtonPoly.h}中再声明一个函数\texttt{Newton\_dif\_at\_x}来计算导数值.可是在编程里怎么算比较好呢？
	
	注意到我们在编写\texttt{Newton\_poly\_at\_x}时将所求多项式拆分成了如下形式：
	$$ 
		p_n(x)=\bigg(\cdots\Big(a_n\left(x+b_{n-1}\right)+a_{n-1}\Big)(x+b_{n-2})+\cdots+a_1\bigg)(x+b_0)+a_0
	$$
	对其求导，理论上也具有一定的递推性.我们记
	$$ 
		y_{i}(x):=\bigg(\cdots\Big(a_n\left(x+b_{n-1}\right)+a_{n-1}\Big)(x+b_{n-2})+\cdots+a_{i+1}\bigg)(x+b_i)+a_i
	$$
	则有
	$$
		y_n(x)=a_n,y_{i}(x)=y_{i+1}(x)(x+b_i)+a_i,i=0,1,\cdots,n-1.(y_0(x)=p_n(x))
	$$
	根据求导法则，
	\begin{align}
		p_n'(x)&=y_0'(x)=y_1'(x)(x+b_0)+y_1(x)\nonumber\\
		&=\Big(y_2'(x)(x+b_1)+y_2(x)\Big)(x+b_0)+y_1(x)\nonumber\\
		&=\cdots\text{（利用}y_i\text{的递推公式求导，代入导数项直到$y_n'(x)=0$）}\nonumber\\
		&=\bigg(\cdots\Big(y_n(x)\left(x+b_{n-2}\right)+y_{n-1}(x)\Big)(x+b_{n-3})+\cdots+y_2(x)\bigg)(x+b_0)+y_1(x)\nonumber
	\end{align}
	这跟$p_n(x)$具有类似的结构，因此编程就变得容易了.
	
	现在我们来解答这道题.根据上述方法传入数据并得到如下输出结果：
	
	\begin{figure}[H]
		\centering
		\includegraphics[width=0.9\textwidth]{imgs/AnswerD.png}
		\caption{AnswerD}
		\label{D}
	\end{figure}

	这里稍微解释一下“File ProblemD\_data/file.txt created and written successfully.”的意思.在需要描点数据时，我们将输出内容放入\texttt{txt}文件中，而这句话表示数据已成功放入\texttt{txt}文件中.如果没有成功放入，则会显示“Error creating file: ProblemD\_data/file.txt”.这些文件方面的显示对题B、C、E均适用.
	
	回到原题.我们在图\ref{D}中已经计算出了车在$t=10s$时的位置和速度，接下来需要判断车子是否超速.为了判断之，我们在区间$[0,13]$上进行“描点连线”，作出图像$v=v(t)$.将它跟$v=81$图像放在一起，如下图所示.可以很明显地看出，在$t=6s$后的一段时间有轻微超速，而在$t=12s$前后更是严重超速.当然这只是模拟，真实情况可能不会这样（就像题B的Runge现象一样），存在误差是正常的.

	\begin{tikzpicture}  
		\begin{axis}[  
			title = 题D：Hermite插值法模拟车速,
			legend pos=outer north east,
			xlabel = $t$,
			ylabel = {$v$},
			xtick={0, 2, 4, 6, 8, 10, 12, 14},
			ytick={40, 60, 80, 100, 120},
			]  
			\addplot[blue, mark=none] table {ProblemD_data/file.txt};  
			\addlegendentry{$v(t)$}
			\addplot[red, thick] coordinates {(0, 81) (13, 81)};
			\addlegendentry{$v=81$}
		\end{axis}  
	\end{tikzpicture}
	
	\section{ProblemE}
	
	
	\begin{sloppypar}
		(a)这里相当于给出了插值点（天数，记为$d$）$d_0=0,d_1=6,d_2=10,d_3=13,d_4=17,d_5=20,d_6=28$以及两组对应函数值（幼虫平均重量，记为$w=w(d)$）
		$w_1=\{w_{10},\cdots, w_{16}\}=\{6.67,\ 17.3,\ 42.7,\ 37.3,\ 30.1,\ 29.3,\ 28.7\},\ w_2=\{w_{20},\cdots,w_{26}\}=\{6.67,\ 16.1,\ 18.9,\ 15.0,\ 10.6,\ 9.44,\ 8.89\}.$对每个函数值组，结合插值点，由Newton插值法绘制出插值多项式的图像.只需使用\texttt{NewtonPoly.h}中的\texttt{Newton\_poly\_at\_x}，传入插值点，函数值组以及采样点$x$，即可进行描点连线.本题在区间$[0,28]$上使用了281个采样点，最终绘出以下的两条曲线.
		
	\end{sloppypar}
	
	
	
	\begin{tikzpicture}  
		\begin{axis}[  
			title = 题E(a)：绘制每个样本的平均重量曲线,
			legend pos=outer north east,
			xlabel = {Day},
			ylabel = {Average Weight},
			ymajorgrids = true,
			grid style = dashed,
			]  
			\addplot[blue, mark=none] table {ProblemE_data/file1.txt};  
			\addlegendentry{Sp1}
			\addplot[red, mark=none] table {ProblemE_data/file2.txt};
			\addlegendentry{Sp2}
		\end{axis}  
	\end{tikzpicture}

	(b)对于预测15天后样本的幼虫是否死亡，只需要知道其平均质量是否会小于0.根据表格给出的数据，直观上我们感知，幼虫平均重量理论上是在第28天后越来越低，最终变为0的.
	
	我们仿照(a)题，将区间延长到$[0,43]$，使用431个采样点，绘制出以下的两条曲线.
	
	\begin{tikzpicture}  
		\begin{axis}[  
			title = 题E(b)：预测后15天的平均重量曲线,
			legend pos=outer north east,
			xlabel = {Day},
			ylabel = {Average Weight},
			ymajorgrids = true,
			grid style = dashed,
			]  
			\addplot[blue, mark=none] table {ProblemE_data/file1_predict.txt};  
			\addlegendentry{Sp1}
			\addplot[red, mark=none] table {ProblemE_data/file2_predict.txt};
			\addlegendentry{Sp2}
		\end{axis}  
	\end{tikzpicture}
	
	可以看到，在第28天过后，两组的幼虫平均重量均显著上升，样本1的重量在第43天更是达到了$10^4$的量级，这很明显是不符合直观乃至实际现象的，主要原因是插值点给的信息太少（或者是密度分布不合理），导致插值多项式给出的误差较大.另外，我们根据先前提到的思路写出代码，运行结果如图\ref{E}.
	
	\begin{figure}[H]
		\centering
		\includegraphics[width=0.9\textwidth]{imgs/AnswerE.png}
		\caption{AnswerE}
		\label{E}
	\end{figure}

	我们看到，输出显示“在第0天时，样本1的幼虫已死亡”.这是因为样本1的插值多项式误差太大的缘故，在前5天基本上重量是一个负值.而样本2也可能只是在区间$[0,28]$上拟合较好，后续也有很大的误差（如输出结果，达到了$10^3$的数量级，不符实际）.
	
\end{document}